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Abstract 

o 

The most widely applied strategy for workload sharing is to equalize the workload assigned to each resource. In 
^ mobile multi-agent systems, this principle directly leads to equitable partitioning policies in which (i) the workspace 

is divided into subregions of equal measure, (ii) there is a bijective correspondence between agents and subregions, 
and (iii) each agent is responsible for service requests originating within its own subregion. In this paper, we design 
CO provably correct, spatially-distributed and adaptive policies that allow a team of agents to achieve a convex and 

equitable partition of a convex workspace, where each subregion has the same measure. We also consider the issue 
of achieving convex and equitable partitions where subregions have shapes similar to those of regular polygons. Our 
approach is related to the classic Lloyd algorithm, and exploits the unique features of power diagrams. We discuss 
^ possible applications to routing of vehicles in stochastic and dynamic environments. Simulation results are presented 

I— I and discussed. 

> 

I. Introduction 

(S| In the near future, large groups of autonomous agents will be used to perform complex tasks including transporta- 

tion and distribution, logistics, surveillance, search and rescue operations, humanitarian demining, environmental 

cn 

O monitoring, and planetary exploration. The potential advantages of multi-agent systems are, in fact, numerous. For 

On 

Q instance, the intrinsic parallelism of a multi-agent system provides robustness to failures of single agents, and in 
• • 

^ many cases can guarantee better time efficiency. Moreover, it is possible to reduce the total implementation and 

^ Operation cost, increase reactivity and system reliability, and add flexibility and modularity to monolithic approaches. 

H 

In essence, agents can be interpreted as resources to be shared among customers. In surveillance and exploration 
missions, customers are points of interests to be visited; in transportation and distribution applications, customers 
are people demanding some service (e.g., utility repair) or goods; in logistics tasks, customers could be troops in 
the battlefield. Finally, consider a possible architecture for networks of autonomous agents performing distributed 
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sensing: a set of n cheap sensing devices (sensing nodes), distributed in the environment, provides sensor mea- 
surements, while m sophisticated agents (cluster heads) collect information from the sensing nodes and transmit it 
(possibly after some computation) to the outside world. In this case, the sensing nodes represent customers, while 
the agents, acting as cluster heads, represent resources to be allocated. 

The most widely applied strategy for workload sharing among resources is to equalize the total workload assigned 
to each resource. While, in principle, several strategies are able to guarantee workload-balancing in multi-agent 
systems, equitable partitioning policies are predominant Q, (21, O, (H. A partitioning policy is an algorithm that, 
as a function of the number m of agents and, possibly, of their position and other information, partitions a bounded 
workspace A C into m openly disjoint regions A^, for z G {1, . . . , m}. (Voronoi diagrams are an example of a 
partitioning policy.) In the resource allocation problem, each agent i is assigned to subregion A^, and each customer 
in Ai receives service by the agent assigned to Ai. Accordingly, if we model the workload for subregion S ^ A 
as \s = X{x) dx, where X{x) is a measure over A, then the workload for agent i is A^^. Given this preface, 
load-balancing calls for equalizing the workload A^^ in the m subregions or, in equivalent words, to compute an 
equitable partition of the workspace A (i.e., a partition where A^^ = A^/^, for all 0- 

Equitable partitioning policies are predominant for three main reasons: (i) efficiency, (ii) ease of design and (iii) 
ease of analysis. Equitable partitioning policies are, therefore, ubiquitous in multi-agent system applications. To 
date, nevertheless, to the best of our knowledge, all equitable partitioning policies inherently assume a centralized 
computation of the workspace partition. This fact is in sharp contrast with the desire of a fully distributed architecture 
for a multi-agent system. The lack of a fully distributed architecture limits the applicability of equitable partitioning 
policies to limited-size multi-agent systems operating in a known static environment. 

The contribution of this paper is three-fold. First, we design provably correct, spatially-distributed, and adaptive 
policies that allow a team of agents to achieve a convex and equitable partition of a convex workspace. Our 
approach is related to the classic Lloyd algorithm from vector quantization theory 0, 0, and exploits the unique 
features of power diagrams, a generalization of Voronoi diagrams (see |7 | for another interesting application of 
power diagrams in mobile sensor networks). Second, we provide extensions of our algorithms to take into account 
secondary objectives, as for example, control on the shapes of the subregions. Our motivation, here, is that equitable 
partitions in which subregions are thin slices are, in most applications, impractical: in the case of dynamic vehicle 
routing, for example, a thin slice partition would directly lead to an increase in fuel consumption. Third, we 
discuss some applications of our algorithms; we focus, in particular, on the Dynamic Traveling Repairman Problem 
(DTRP) 1 1|, where equitable partitioning policies are indeed optimal under some assumptions. 

Finally, we mention that our algorithms, although motivated in the context of multi-agent systems, are a novel 
contribution to the field of computational geometry. In particular we address, using a dynamical system framework, 
the well-studied equitable convex partition problem (see lO and references therein); moreover, our results provide 
new insights in the geometry of Voronoi diagrams and power diagrams. 

The paper is organized as follows. In Section II we provide the necessary tools from calculus, degree theory and 
geometry. Section III contains the problem formulation, while in Section IV we present preliminary algorithms for 
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equitable partitioning based on leader-election, and we discuss their limitations. Section V is the core of the paper: 
we first prove some existence results for power diagrams, and then we design provably correct, spatially-distributed, 
and adaptive equitable partitioning policies that do not require any leader election. In Section VI we extend the 
algorithms developed in Section V to take into account secondary objectives. Section VII contains simulations 
results. Finally, in Section VIII, we provide an application of our algorithms to the DTRP problem, and we draw 
our conclusions. 

II. Background 

In this section, we introduce some notation and briefly review some concepts from calculus, degree theory and 
geometry, on which we will rely extensively later in the paper. 

A. Notation 

Let II • II denote the Euclidean norm. Let A be a compact, convex subset of R^. We denote the boundary of A as 
dA and the Lebesgue measure of A as \A\. We define the diameter of A as: diameter (A) = sup{| |p— g| | | g G A}. 
The distance from a point x to a set M is defined as dist(x, M) = inf^^M ||^ — p||- We define = {1, 2, • • • , m}. 
Let G = {gij - ■ ■ jgm) C denote the location of m points. A partition (or tessellation) of A is a collection of 
m closed subsets A = {Ai, • • • , Am} with disjoint interiors whose union is A. A partition A = {Ai, • • • , Am} is 
convex if each Ai, i e ImA^ convex. 

Given a vector space V, let F(V) be the collection of finite subsets of V. Accordingly, F(R^) is the collection 
of finite point sets in R^. Let G(R'^) be the set of undirected graphs whose vertex set is an element of F(R'^) (we 
assume the reader is familiar with the standard notions of graph theory as defined, for instance, in |9, Chapter 1]). 

Finally, we define the saturation function sata,6(x), with a < b, as: 

1 if x>b 

ssita,b{x) = ^ {x-a)/{b-a) if a < x < b (1) 
otherwise 

B. Variation of an Integral Function due to a Domain Change. 

The following result is related to classic divergence theorems ifTOl . Let Q = Q{y) C A be a region that depends 
smoothly on a real parameter y e R and that has a well-defined boundary dQ{y) for all y. Let be a density 
function over A. Then 

h{x) dx = I - n{x)\ h{x) dx ^ (2) 

JdQ{y)^dy ^ 

where v • w denotes the scalar product between vectors v and w, where n{x) is the unit outward normal to dVt{y), 
and where dx/dy denotes the derivative of the boundary points with respect to y. 
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C. A Basic Result in Degree Theory 

In this section, we state some results in degree theory that will be useful in the remainder of the paper. For a 
thoroughly introduction to the theory of degree we refer the reader to ifTTI . 

Let us just recall the simplest definition of degree of a map /. Let / : X ^ F be a smooth map between 
connected compact manifolds X and Y of the same dimension, and let p G F a regular value for / (regular values 
abound due to Sard's lemma). Since X is compact, f~^{p) = {xi, . . . , x^} is a finite set of points and since p is a 
regular value, it means that /f/. : Ui f{Ui) is a local diffeomorphism, where /7i is a suitable open neighborhood 
of Xi. Diffeomorphisms can be either orientation preserving or orientation reversing. Let (i+ be the number of points 
Xi in f~^{p) at which / is orientation preserving (i.e. det(Jac(/)) > 0, where Jac(/) is the Jacobian matrix of 
/) and d~ be the number of points in f~^{p) at which f is orientation reversing (i.e. det(Jac(/)) < 0). Since X 
is connected, it can be proved that the number (i+ — d~ is independent on the choice ofpeY and one defines 
deg(/) := d'^ — d~ . The degree can be also defined for a continuous map f : X ^ Y among connected oriented 
topological manifolds of the same dimensions, this time using homology groups or the local homology groups at 
each point in f~^{p) whenever the set f~^{p) is finite. For more details see ifTTIl . 

The following result will be fundamental to prove some existence theorems and it is a direct consequence of the 
theory of degree of continuous maps among spheres. 

Theorem 2.1: Let / : be a continuous map from a closed n-ball to itself. Call S^~^ the boundary of 

B^, namely the (n — l)-sphere and assume that fs^^ : S'^ ^ S'^ di map with deg(/) ^ 0. Then / is onto B^ . 

Proof: Since / as a map from S^~^ to S^~^ is different from zero, then the map fs^-^ is onto the sphere. 
If / is not onto B^, then it is homotopic to a map B^ S^~^, and then /^r^-i : S^~^ S^~^ is homotopic 
to the trivial map (since it extends to the ball). Therefore /^r^-i : S^~^ S^~^ has zero degree, contrary to the 
assumption that it has degree different from zero. ■ 
In the sequel we will need also the following: 

Lemma 2.2: Let f : ^ a. continuous bijective map from the n-dimensional sphere to itself (n > 1). Then 
deg(/) = ±1. 

Proof: The map / is a continuous bijective map from a compact space to a Hausdorff space, and therefore it 
is a homeomorphism. Now, a homeomorphism f : S'^ ^ S'^ has degree ±1 (see, for instance, (TTl Page 136]). ■ 

D. Voronoi Diagrams and Power Diagrams 

We refer the reader to |[T2ll and ifTSll for comprehensive treatments, respectively, of Voronoi diagrams and power di- 
agrams. Assume, first, that G is an ordered set of distinct points. The Voronoi diagram V(G) = {Vi{G), • • • , Vm{G)) 
of A generated by points (^i, • • • , ^m) is defined by 

V,{G) = {xeA\ \\x-g,\\ < Wx-gjl^jy^ijelm}- (3) 

We refer to G as the set of generators of V(G), and to Vi{G) as the Voronoi cell or region of dominance of the 
i-th generator. For gi^gj G G, i 7^ j, we define the bisector between gi and gj as b{gi^gj) = {x G A\ \\x — gi\\ = 
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||x — The face b{gi^gj) bisects the Hne segment joining gi and gj, and this Hne segment is orthogonal to the 
face (Perpendicular Bisector Property). The bisector divides A into two convex subsets, and leads to the definition 
of the set D{gi^gj) = {x e A\ \\x — gi\\ < \\x — gj\\}; we refer to D{gi^gj) as the dominance region of gi over 
gj. Then, the Voronoi partition V{G) can be equivalently defined as Vi{G) = f]j^j^^^-y D{gi^ gj). This second 
definition clearly shows that each Voronoi cell is a convex set. Indeed, a Voronoi diagram of A is a convex partition 
of A (see Fig. |l(a)| ). The Voronoi diagram of an ordered set of possibly coincident points is not well-defined. We 
define 

Tcoinc = {(^1, - ,gm) ^ A'^\gi= gj for some z / j G {1, • • • , m}}. (4) 

Assume, now, that each point gi ^ G has assigned an individual weight Wi e i e Im'A^t W = {wi^ • • • , Wm)- 
We define the power distance 

dp{x,gi]Wi) = \\x- gif - Wi. (5) 

We refer to the pair (^^, Wi) as a power point. We define Gw = {{dii'^i)^ ' ' ' i {dm^ ^m)^ • Two power points 
(gi^Wi) and (gj^Wj) are coincident if g^ = gj and Wi = Wj. Assume, first, that Gw is an ordered set of distinct 
power points. Similarly as before, the Power diagram V{Gw) = {yi{Gw)r ' ' ^ym{Gw)) of A generated by 
power points (^(^i, ^i), • • • , (^^n, ^m)) is defined by 

V^{Gw) = {xeA\ \\x - g,f -w,< \\x - g^ f - Wj, Vj / j G /^}. (6) 

We refer to Gw as the set of power generators of V{Gw), and to Vi{Gw) as the power cell or region of dominance 
of the i-th power generator; moreover we call g^ and Wi, respectively, the position and the weight of the power 
generator (gi^Wi). Notice that, when all weights are the same, the power diagram of A coincides with the Voronoi 
diagram of A. As before, power diagrams can be defined as intersection of convex sets; thus, a power diagram 
is, as well, a convex partition of A. Indeed, power diagrams are the generalized Voronoi diagrams that have the 
strongest similarities to the original diagrams 11141 . There are some differences, though. First, a power cell might 



be empty. Second, gi might not be in its power cell (see Fig. 1(b)). Finally, the bisector of (gi^Wi) and {gj^Wj), 
i + h is 

= G A| {g.-giYx = \{\g,f - \\g,f ^ Wi - w,)}. (7) 

Hence, b(^{gi^Wi)^{gj^Wj)^ is a face orthogonal to the line segment g^gj and passing through the point g*j given 
by 

ll«.l|2 _ ||„.||2 _i_ _ 



\\9j\\^-\\9if+m-wj 



this last property is crucial in the remaining of the paper: it means that, by changing weights, it is possible to 
arbitrarily move the bisector between the positions of two power generators, while still preserving the orthogonality 
constraint. 

The power diagram of an ordered set of possibly coincident power points is not well-defined. We define 

Tcoinc = I (^{gi.wi), • • • , {gm, Wm)^ G (A X R)"^ I gi = gj and Wi = Wj for some i 7^ j G {1, . . . , m}|. (8) 
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Notice that we used the same symbol as in Eq. (|4]): the meaning will be clear from the context. 

For simplicity, we will refer to Vi{G) (Vi{Gw)) as Vi. When the two Voronoi (power) cells Vi and Vj are 
adjacent (i.e., they share a face), gi ({gi, Wi)) is called a Voronoi (power) neighbor of gj ({gj^Wj)), and vice- versa. 
The set of indices of the Voronoi (power) neighbors of gi {{gi^ wi)) is denoted by Ni. We also define the (z, j)-face 
as /^ij = Vi{^ Vj. 




(a) A Voronoi Diagram. 



Fig. 1. Voronoi diagrams and power diagrams. 



(b) A power diagram. The weights Wi are as- 
sumed positive. Circles represent the magnitudes 
of weights. Power generator {g2,W2) has an 
empty cell. Power generator (g^,ws) is outside 
its region of dominance. 



E. Proximity Graphs and Spatially-Distributed Control Policies for Robotic Networks 

Next, we shall present some relevant concepts on proximity graph functions and spatially-distributed control 
policies; we refer the reader to ifTSl for a more detailed discussion. A proximity graph function Q : F(R^) 
G{R^) associates to a point set V G F(R'^) an undirected graph with vertex set V and edge set Sg{V), where 
Sg : F(R^) ^ F(R'^ x R^) has the property that Sg{V) CV xV\ dmg{V x V) for any V . Here, diag(7^ xV) = 
{{p^p) ^VxV\ p eV}.ln other words, the edge set of a proximity graph depends on the location of its vertices. 
To each proximity graph function, one can associate the set of neighbors map Ng : R^ x F(R^) F(R^), defined 
by 

Ng{p,r) = {qer\ {p,q)e£g{ru{p})}. 

Two examples of proximity graph functions are: 
(i) the Delaunay graph G ^ Gy{G) = {G,Eg^{G)) has edge set 

£qAG) = {{g^.gJ)eGxG\ diag(G x G)\ V,{G) n Vj{G) + 0}, 
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where Vi{G) is the i-th cell in the Voronoi diagram V(G); 
(ii) the power-Delaunay graph Gw ^ Gp{Gw) = {Gw^^g^iGw)) has edge set 



Sg,{Gw) = {(^gi,Wi),{gj,Wj)^ e Gw x Gw\dmg{Gw x Gw)\ Vi{Gw) r\Vj{Gw) + 0}, 



where Vi{Gw) is the i-th cell in the power diagram V(Gvr)- 
We are now in a position to discuss spatially-distributed algorithms for robotic networks in formal terms. Let 
P(t) = . . . ,Pm{t)) G be the ordered set of positions of m agents in a robotic network. We denote the 

state of each agent i G at time t as Xi(t) G ixi{t) can include the position of agent i as well as other 
information). With a slight abuse of notation, let us denote by Ii(t) the information available to agent i at time t. The 
information vector Ii{t) is a subset of x{t) = {xi{t), . . . , Xm{t)) of the form Ii{t) = {xi^ (t), . . . , Xi^{t)}, k < m. 
We assume that Ii{t) always includes Xi{t). Let ^ be a proximity graph function defined over P{t) (respectively 
over Pw{t) if we also consider a weight Wi{t) for each robot i G /m); we define I^^ (t) as the information vector 
with the property Xi{t) G I^^ (t), and, for j ^ i. 



In words, the information vector /• ^ (t) coincides with the states of the neighbors (as induced by Q) of agent i 
together with the state of agent i itself. 

Let = . . . ^fim{Im{t)) be a feedback control policy for the robotic network. The policy /i is 

spatially distributed over Q if for each agent i E Im and for all t 



In other words, through information about its neighbors according to each agent i has sufficient information to 
compute the control /i^. 

III. Problem Formulation 

A total of m identical mobile agents provide service in a compact, convex service region A C R^. Let A be 
a measure whose bounded support is A (in other words, A is not zero only on A); for any set S, we define the 
workload for region S sls Xs = X{x) dx. The measure A models service requests, and can represent, for example, 
the density of customers over A, or, in a stochastic setting, their arrival rate. Given the measure A, a partition {Ai}i 
of the workspace A is equitable if A^. = Xaj for all j E Im- 

A partitioning policy is an algorithm that, as a function of the number m of agents and, possibly, of their position 
and other information, partitions a bounded workspace A into m openly disjoint subregions Ai, i e Im- Then, each 
agent i is assigned to subregion Ai, and each service request in Ai receives service from the agent assigned to 
Ai. We refer to subregion A^ as the region of dominance of agent i. Given a measure A and a partitioning policy, 
m agents are in a convex equipartition configuration with respect to A if the associated partition is equitable and 
convex. 
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In this paper we study the following problem: find a spatially -distributed (in the sense discussed in Section [n]) 
equitable partitioning policy that allows m mobile agents to achieve a convex equipartition configuration (with 
respect to A). Moreover, we consider the issue of convergence to equitable partitions with some special properties, 
e.g., where subregions have shapes similar to those of regular polygons. 

IV. Leader-Election Policies 

We first describe two simple algorithms that provide equitable partitions. A first possibility is to "sweep" A 
from a point in the interior of A using an arbitrary starting ray until A^i = A^/^, continuing the sweep until 
= -^a/^, etc. A second possibility is to slice, in a similar fashion, A. The resulting equitable partitions are 
depicted in Fig. [2] 




(a) Sweeping A (b) Slicing A 



Fig. 2. Equitable partitions by sweeping and slicing (assuming a uniform measure A). 

Then, a possible solution could be to (i) run a distributed leader election algorithm over the graph associated 
to some proximity graph function Q (e.g., the Delaunay graph); (ii) let each agent send its state Xi{t) to the 
leader; (iii) let the leader execute either the sweeping or the slicing algorithms described above; finally, (iv), let 
the leader broadcast subregion's assignments to all other agents. Such conceptually simple solution, however, can 
be impractical in scenarios where the density A changes over time, or agents can fail, since at every parameter's 
change a new time-consuming leader election is needed. Moreover, the sweeping and the slicing algorithms provide 
long and skinny subregions that are not suitable in most applications of interest (e.g., vehicle routing). 

We now present spatially-distributed algorithms, based on the concept of power diagrams, that solve both issues 
at once. 

V. Spatially-Distributed Gradient-Descent Law for Equitable Partitioning 
We start this section with an existence theorem for equitable power diagrams. 



March 30, 2009 



DRAFT 



ffiEE TRANSACTIONS ON AUTOMATIC CONTROL, SUBMITTED FOR PUBLICATION AS A FULL PAPER. 



9 



A. On the Existence of Equitable Power Diagrams 

As shown in the next theorem, an equitable power diagram (with respect to any A) exists for any vector of 
distinct points G = (^i, . . . , in A. 

Theorem 5.1: Let A be a bounded, connected domain in R^, and A be a measure on A. Let G = (^i, . . . ^Qn) 
be the positions of n > 1 distinct points in A. Then, there exist weights Wi, i G /n, such that the power points 
(^(^1, ^1 ),•••, (^n, ^n)^ generate a power diagram that is equitable with respect to A. Moreover, given a vector of 
weights w* that yields an equitable power diagram, the set of all vectors of weights yielding an equitable power 
diagram is = {^* + . . . , 1]}, with t G M. 

Proof: It is not restrictive to assume that A^ = 1 (i.e., we normalize the measure of A), since A is bounded. 
The strategy of the proof is to use a topological argument to force existence. 

First, we construct a weight space. Let D = diameter(74), and consider the cube C := [—D^D]'^. This is the 
weight space and we consider weight vectors W taking value in C. Second, consider the standard n-simplex of 
measures A^i , • • • , A^^ (where Ai, . . . , are, as usual, the regions of dominance). This can be realized in as 
the subset of defined by Y^^^i A^. = 1 with the condition A^. > 0. Let us call this set "the measure simplex A' 
(notice that it is (n — 1) -dimensional). 

There is a map f : C ^ A associating, according to the power distance, a weight vector W with the corresponding 
vector of measures (A^i, • • • , A^^). Since the points in G are assumed to be distinct, this map is continuous. 

We will now use induction on n, starting with the base case n = 3 (the statement for n = 1 and n = 2 is trivially 
checked). We study in detail the case for n = 3, where visualization is easier, in order to make the proof more 
transparent. When n = 3, the weight space C is a three dimensional cube with vertices = [—D^ —D^ —D]^ vi = 
[D,-D, -D] ,V2 = [-D,D, -D] ,vs = [-D,-D,D],v^ = [D,-D,D],v^ = [-D,D,D],VG = [D,D, -D] and 
vy = [D^D^D]. The measure simplex A is, instead, a triangle with vertices Ui^U2j us that correspond to the cases 
A^i = 1, Aa2 = 0, A^s = 0, A^i = 0, = 1, = 0, and A^i = 0, = 0, ^As = 1, respectively. Moreover, 
call 61,62 and 63 the edges opposite the vertices ui^U2^us respectively. The edges 6^ are, therefore, given by the 
condition A^. G 6^ <^ A^. = 0. 

Let us return to the map f : C ^ A (now in the case of three generators). Observe that the map / sends vq 
the unique point po of A corresponding to the measures of usual Voronoi cells (since the weights are all equal). 
Call li the edge ^o^i; then, it is immediate to see that the image of h through / is a path 71 in A joining po to 
ui. Analogously, the image of I2 = V0V2 through / is a path 72 in A joining po to U2 and, finally, the image of 
^3 = ^0^3 through / is a path 73 connecting po to us (see Fig. |3]). Now, we observe that paths {7i}{i=i,2,3} do 
not intersect except in pq. To prove this, start by observing that the image through / of all the points on the main 
diagonal of the cube joining vq with V7 is equal to a single point po G A. This is due to the fact that only the 
differences among weights change the vector (A^^ , A^2 5 ^^3), i-^., if all weights are increased by the same quantity, 
the vector (A^i , A^s , ^^3) does not change. We will prove this in detail for the case of n-generators in the next few 
paragraphs. In particular, the image of the diagonal vqVy is exactly the point for which the measures are those of a 
Voronoi partition. Now let us understand what are the "fibers" of /, that is to say, the loci where / is constant. Since 
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the measures of the regions of dominance do not change if the differences among the weights are kept constant, 
then the fibers of / in the weight space C are given by the equations wi—W2 = ci and W2 — ws = C2. Rearranging 
these equations, it is immediate to see that wi = ws -\- ci -\- C2, W2 = ws -\- C2 and ws = ws, therefore taking ws 
as parameter we see that the fibers of / are straight fines parallel to the main diagonal vqVj. On the weight space 
C let us define the following equivalence relation: w = w' if and only if they are on a line parallel to the main 
diagonal vqVj. Map f : C ^ A induces a continuous map (still called / by abuse of notation) from C/ = to ^ 
having the same image. Let us identify C/ =. It is easy to see that any line in the cube parallel to the main diagonal 
vqVj is entirely determined by its intersection with the three faces F3 = {w^, = —D} DC, F2 = {w2 = —D} DC 
and Fi = {wi = —D} n C. Call JT the union of these faces. We therefore have a continuous map f : ^ A that 
has the same image of original /; besides, in the few next paragraphs we will prove in general (i.e for the case 
with n generators) that the induced map f : !F ^ A is injective by construction. 

Observe that JT is homeomorphic to 5^, the 2-dimensional ball, like A itself. Up to homeomorphisms, therefore, 
the map f : ^ A can be viewed as a map (again called / by abuse of notation), f : B'^ ^ B'^. Consider 
the closed loop a given by ^2^5, V5V3, v^vi, viVe, ^6^2 with this orientation. This loop is the boundary of 

JT and we think of it also as the boundary of B'^. It is easy to see that / maps a onto the boundary of A, and 
since / is injective, the inverse image of any point in the boundary of A is just one element of a. Identifying the 
boundary of A with (up to homeomorphisms) and the loop a with (up to homeomorphisms) we have a 
bijective continuous map /51 : ^ 5*^. By Lemma ( |2.2| ) deg(/) = ±1 and therefore / is onto A, using Theorem 



Fig. 3. Construction used for the proof of existence of equitable power diagrams. 

Now we extend the same idea to the case of n generators and we will use also induction on the number of 
agents. Therefore, we suppose that we have proved that the map / is surjective for n — 1 agents and we show how 
to use this to show that the map is surjective for n agents. 



Vq = (-D,-D,-D) 
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If we have n generators, the weight space is given by an n dimensional cube C = [—D^D]^, in complete analogy 
with the case of 3 generators. The n simplex of the areas A is again defined as a {(A^i, • • • , A^^) G M^} such 
that Xa^ > for z = 1, . . . , n and Yl7=i -^^^ ~ ^' Notice that A is homeomorphic to the (n — 1) -dimensional 
ball B^~^. As before we have a continuous map / : C ^ ^. It is easy to see that / is constant on the sets of 
the form := {{w* + t(l,...,l)} flC, t G M}, that is whenever two sets of weights differ by a common 
quantity, they are mapped to the same point in A. Moreover, fixing a point Q G ^ we have that f~^{Q) is given 
just by a set of the form for a suitable w*. Indeed, assume this is not the case, then the vector of measures 
(A^^, . . . , A^^) is obtained via / using two sets of weights: := {wl^ . . . w}^) and := {wf^ . . . w'^), and 
and don't belong to the same W^, namely it is is not possible to obtain as + . . . , 1) for a 
suitable t. This means that the vector difference — is not a multiple of (1, . . . , 1). Therefore, there exists a 
nonempty set of indexes J, such that w'j —Wj >w1 — w\, whenever j G J and for all A: G {1, ... n} and such that 
the previous inequality is strict for at least one /c*. Now among the indexes in J, there exists at least one of them, 
call it j* such that the agent j* is a neighbor of agent A:*, due to the fact that the domain A is connected. But since 
— w^j^ > wl^ — wl^, and Wj^ — Wj^ > wl — wl for all A: G {1, ... , n}, this implies that the measure A^^* 
corresponding to the choice of weights is strictly greater that Xaj* corresponding to the choice of weights W^. 
This proves that f~^{Q) is given only by sets of the form W^. 

We introduce an equivalence relation on C, declaring that two sets of weights and are equivalent if 
and only if they belong to the same W^. Let us call = this equivalence relation. It is immediate to see that / 
descends to a map, still called / by abuse of notation, f : C/ =^ A and that / is now injective. It is easy also 
to identify C/ = with the union of the (n — 1) -dimensional faces of C given by = U^^i(C H {wi = —D}). In 
this way we get a continuous injective map f : ^ A that has the same image as the original /. Notice also that 
JT is homeomorphic to the closed (n — 1) -dimensional ball, so up to homeomorphism / can be viewed as a map 

We want to prove that the map /^jr, given by the restriction of / to 9jF is onto dA. To see this, consider one 
of the (n — 2) -dimensional faces dAi of dA, which is identified by the condition A^^ = 0. Consider the face Fi 
in JT, where Fi is given by Fi := C D {wi = —D}. We claim that the Si := dFi H dJ^ is mapped onto dAi by /. 
Observe that the Si is described by the following equations Si = Uj^i{{wi = —D^Wj = D}nJ^), so Si is exactly 
equivalent to a set of type for n — 1 agents. Moreover observe that dAi can also be identified with the measure 
simplex for n — 1 agents. By inductive hypothesis therefore, the map f : Si ^ dAi is surjective, and therefore 
also the map foj^ is onto dA. Since /ajr is a bijective continuos map among (n — 2) -dimensional spheres, (up to 
homeomorphism), it has degree ±1 by Lemma ( |2.2| ). Finally we conclude that / is onto A, using again Theorem 

■ 

Some remarks are in order. 

Remark 5.2: The above theorem holds for any bounded, connected domain in R^. Thus, the case of a compact, 
convex subset of is included as a special case. Moreover, it holds for any measure A absolutely continuous with 
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respect to the Lebesgue measure, and for any vector of distinct points in A. 

Remark 5.3: In the proof of the above theorem, we actually proved that for any measure vector {XAi}i=i,...m in 
A, there exists a weight vector w e C realizing it through the map /. This could be useful in some applications. 

Remark 5.4: Since all vectors of weights in W yield exactly the same power diagram, we conclude that the 
positions of the generators uniquely induces an equitable power diagram. 



B. State, Region of Dominance, and Locational Optimization 

The first step is to define the state for each agent i. We let Xi{t) be the power generator {gi{t)^ Wi{t)) e A x R, 
where gi{t) = Pi{t) (i.e., the position of the power generator coincides with the position of the agent) We, then, 
define the region of dominance for agent i as the power cell Vi = Vi{Gw), where Gw = (^(^i, ^i), • • • , (^m, ^m)^ • 
We refer to the partition into regions of dominance induced by the set of generator^ as V{Gw)- 

In light of Theorem |5.1[ the key idea is to enable the weights of the generators to follow a spatially-distributed 
gradient descent law (while maintaining the positions of the generators fixed) such that an equitable partition is 
reached. Assume, henceforth, that the positions of the generators are distinct, i.e., gi 7^ gj for i j. Define the set 

S = {{wi,...,Wm)eR^\Xv, >0 \/ieIm}- (9) 

Set S contains all possible vectors of weights such that no region of dominance has measure equal to zero. 
We introduce the locational optimization function Hy • S ^ M>o^ 

where W = {wi^ • • • , Wm)- 

Lemma 5.5: A vector of weights that yields an equitable power diagram is a global minimum of Hy- 
Proof: Consider the relaxation of our minimization problem: 



mm 

3^1 , • • • ,Xm 



^ ; s.t. = a > 0, Xi > 0, i e Im, 



=1 



where the linear equality constraint follows from the fact that X^I^i Iv {W) ^{^)^^ = Sa ^(^) = ^ ^^d where 
the constant a is greater than zero since A is a measure whose bounded support is A. By Lagrange multiplier 
arguments, it is immediate to show that the global minimum for the relaxation is Xi = a/m for all i. Since 



Theorem 5.1 establishes that there exists a vector of weights in S that yields an equitable power diagram, we 



conclude that this vector is a global minimum of Hy. 

^Henceforth, we assume that A is a compact, convex subset of R-^. 
^For short, we will refer to a power generator simply as a generator. 
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C. Smoothness and Gradient of Hy 

We now analyze the smoothness properties of Hy. In the following, let = \\gj — gi\\. 

Theorem 5.6: Assume that the positions of the generators are distinct, i.e., gi 7^ gj for i ^ j. Given a measure 
A, the function Hy is continuously differentiable on S, where for each i G {1, . . . , m} 

Furthermore, the critical points of Hy are vectors of weights that yield an equitable power diagram. 

Proof: By assumption, g^ ^ gj for i / j, thus the power diagram is well defined. Since the motion of a weight 
Wi affects only power cell Vi and its neighboring cells Vj for j e Ni, we have 

dHy ^ _^dXy^ _ y ^dXy^ ^^^^ 

Now, the result in Eq. ([2]) provides the means to analyze the variation of an integral function due to a domain 
change. Since the boundary of Vi satisfies dVi = UjA^j U Bi, where Aij = Aji is the edge between Vi and Vj, 
and Bi is the boundary between Vi and A (if any, otherwise Bi = 9), we have 

=0 

where we defined riij as the unit normal to Aij outward of Vi (therefore we have riji = —riij). The second term is 
trivially equal to zero if 5^ = 0; it is also equal to zero if Bi ^ 0, since the integrand is zero almost everywhere. 
Similarly, 

SAv. [ f dx 



^ • nji{x) \ \{x) dx. (14) 



dwi J /\.\ dwi 

To evaluate the scalar product between the boundary points and the unit outward normal to the border in Eq. 
(13) and in Eq. ([14]), we differentiate Eq. ([T]) with respect to Wi at every point x G Aij\ we get 

^^^^ 

From Eq. ([T]) we have Uij — {gj — gi) /Wdj — 9i\\, and the desired explicit expressions for the scalar products in 
Eq. ( p3] ) and in Eq. ([14]) follow immediately (recalling that riji = —riij). 

Collecting the above results, we get the partial derivative with respect to Wi. The proof of the characterization of 
the critical points is an immediate consequence of the expression for the gradient of Hy; we omit it in the interest 
of brevity. ■ 



Remark 5.7: The computation of the gradient in Theorem 5.6 is spatially-distributed over the power-Delaunay 
graph, since it depends only on the location of the other agents with contiguous power cells. 

Example 5.8 (Gradient of Hy for uniform measure): The gradient of Hy simplifies considerably when A is 

constant. In such case, it is straightforward to verify that (assuming that A is normalized) 

dHy 



A\ 7,AIK-P WA^J' 



9 m 2|A|^^^7,, 



(16) 



where 5ij is the length of the border 



March 30, 2009 



DRAFT 



ffiEE TRANSACTIONS ON AUTOMATIC CONTROL, SUBMITTED FOR PUBLICATION AS A FULL PAPER. 



14 



D. Spatially -Distributed Algorithm for Equitable Partitioning 

Consider the set U = . . . , Wm) ^ \ Xli^i = o|- Since adding an identical value to every weight 

leaves all power cells unchanged, there is no loss of generality in restricting the weights to U; let Q = S f) U. 
Assume the generators' weights obey a first order dynamical behavior described by 



Consider Hy an objective function to be minimized and impose that the weight Wi follows the gradient descent 
given by ( pT| ). In more precise terms, we set up the following control law defined over the set Q 

where we assume that the partition V(iy) = {^i, . • . , Kn} is continuously updated. One can prove the following 
result. 

Theorem 5.9: Assume that the positions of the generators are distinct, i.e. gi 7^ gj for i j. Consider the 
gradient vector field on ft defined by equation ^Tf) . Then generators' weights starting at t = at W {0) e ft, and 
evolving under ([17]) remain in ft and converge asymptotically to a critical point of Hy, i.e., to a vector of weights 
yielding an equitable power diagram. 

Proof: We first prove that generators' weights evolving under (Tf) remain in ft and converge asymptotically 
to the set of critical points of Hy- By assumption, g^ ^ gj for i j, thus the power diagram is well defined. First, 
we prove that set ft is positively invariant with respect to (Vf) . Recall that ft = S nU. Noticing that control law 
^rf\ is a gradient descent law, we have 

K!(W(t)) < Hv{W{t)) < Hv{W{0)), i e t > 0. 

Since the measures of the power cells depend continuously on the weights, we conclude that the measures of all 
power cells will be bounded away from zero; thus, the weights will belong to S for all t > 0, that is, W{t) G S 
Vt > 0. Moreover, the sum of the weights is invariant under control law (Vf) . Indeed, 

o ^ ^TT ^ 



dt ^ dwi ^ ^ 2jij VA?/ a?/ . 

i=i * i=i jeNi '^^ 

since jij = jji, Aij = Aji, and j e Ni ^ i e Nj. Thus, we have W{t) eU \/t>0. Since W{t) G 5 Vt > and 
W{t) G Vt > 0, we conclude that W{t) e S DU = ft, Wt > 0, that is, set ft is positively invariant. 

Second, Hy ' ft M>o is clearly non-increasing along system trajectories, that is, Hy < in ft. 

Third, all trajectories with initial conditions in ft are bounded. Indeed, we have already shown that XlHi = 
along system trajectories. This implies that weights remain within a bounded set: If, by contradiction, a weight 
could become arbitrarily positive large, another weight would become arbitrarily negative large (since the sum of 
weights is constant), and the measure of at least one power cell would vanish, which contradicts the fact that S is 
positively invariant. 
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Finally, by Theorem |5.6[ Hy is continuously differentiable in Q.. Hence, by invoking the LaSalle invariance 



principle (see, for instance, |6|), under the descent flow ([17]), weights will converge asymptotically to the set of 



critical points of Hy, that is not empty as confirmed by Theorem 5.1 



Indeed, by Theorem 5.1 we know that all vectors of weights yielding an equitable power diagram differ by 
a common translation. Thus, the largest invariant set of Hy in Vt contains only one point. This implies that 
lim^^oo W{t) exists and it is equal to a vector of weights that yields an equitable power diagram. 

■ 

Some remarks are in order. 



Remark 5.10: By Theorem 5.9 for any set of generators' distinct positions, convergence to an equitable power 
diagram is global with respect to Vt. Indeed, there is a very natural choice for the initial values of the weights. 
Assuming that at t = agents are in A and in distinct positions, each agent initializes its weight to zero. Then, 
the initial partition is a Voronoi tessellation; since A is positive on A, each initial cell has nonzero measure, and 
therefore 1^(0) G ^ (the sum of the initial weights is clearly zero). 

Remark 5.11: The partial derivative of Hy with respect to the i-\h weight only depends on the weights of 
the agents with neighboring power cells. Therefore, the gradient descent law ([TT]) is indeed a spatially -distributed 
control law over the power-Delaunay graph. We mention that, in a power diagram, each power generator has an 
average number of neighbors less than or equal to six |[T4ll : therefore, the computation of gradient ([17]) is scalable 
with the number of agents. 

Remark 5.12: The focus of this paper is on equitable partitions. Notice, however, that it is easy to extend the 
previous algorithm to obtain a spatially-distributed (again over the power-Delaunay graph) control law that provides 
any desired measure vector {A^^}. In particular, assume that we desire a partition such that 

where (3i G (0, 1), Xl^i (3j = 1. If we redefine Hy : S ^ ]R>o as 



HviW) = J2 



i=l ^V.(W) ' 



then, following the same steps as before, it is possible to show that, under control law 



Wi 



the solution converges to a vector of weights that yields a power diagram with the property A^^ = A'^a (whose 
existence is guaranteed by Remark [53] ). 

Remark 5.13: Define the set T = A^ \ Fcoinc (where Fcoinc is defined in Eq. ([4])). It is of interest to define and 
characterize the vector- valued function W* : T ft that associates to each non-degenerate vector of generators' 
positions a set of weights such that the corresponding power diagram is equitable. Precisely, we define VK*(G) as 
= lim^^oo W{t)^ where W{t) = (i^i(t), . . . , Wm{t)) is the solution of the differential equation, with fixed 
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vector of generators' positions G, 



Wi(t) = -^{W{t)), wi{0) = 0, ie Im^ 

OWi 



By Theorem 5.9 VK*(G) is a well-defined quantity (since the limit exists) and corresponds to an equitable power 
diagram. The next lemma characterizes W*{G). 
Lemma 5.14: The function W* is continuous on F. 

Proof: See Appendix. ■ 



E. On the Use of Power Diagrams 

A natural question that arises is whether a similar result can be obtained by using Voronoi diagrams (of which 
power diagrams are a generalization). The answer is positive if we constrain A to be constant over A, but it is 
negative for general measures A, as we briefly discuss next. 

Indeed, when A is constant over A, an equitable Voronoi diagram always exists. We prove this result in a slightly 
more general set-up. 

Definition 5.15 (Unimodal Property): Let A C be a bounded, measurable set (not necessarily convex). We 
say that A enjoys the Unimodal Property if there exists a unit vector v G such that the following holds. For 
each 5 G M, define the slice A^ = {x e A^ v ■ x = s}^ and call V^(s) = md-i{A^) the {d — 1) -dimensional 
Lebesgue measure of the slice. Then, the function is unimodal. In other words, attains its global maximum at 
a point 5, is increasing on (— oo, 5], and decreasing on [5, 00). 

The Unimodal Property (notice that every compact, convex set enjoys such property) turns out to be a sufficient 
condition for the existence of equitable Voronoi diagrams for bounded measurable sets (with respect to constant A). 

Theorem 5.16: If A cR^ is any bounded measurable set satisfying the Unimodal Property and A is constant on 
A, then for every m > 1 there exists an equitable Voronoi diagram with m (Voronoi) generators. 

Proof: See Appendix. ■ 

Then, an equitable Voronoi diagram can be achieved by using a gradient descent law conceptually similar to the 
one discussed previously (the details are presented in 11161 ). We emphasize that the existence result on equitable 
Voronoi diagrams seems to be new in the rich literature on Voronoi tessellations. 

While an equitable Voronoi diagram always exists when A is constant over A, in general, for non-constant A, an 
equitable Voronoi diagram fails to exist, as the following counterexample shows. 

Example 5.17 (Existence problem on a line): Consider a one-dimensional Voronoi diagram. In this case a Voronoi 
cell is a half line or a line segment (called a Voronoi line), and Voronoi vertices are end points of Voronoi lines. 
It is easy to notice that the boundary point between two adjacent Voronoi lines is the mid-point of the generators 
of those Voronoi lines. Consider the measure A in Fig. |4j whose support is the interval [0, 1]. Assume m = 5. Let 
6i (z = 1 , . . . , 4) be the position of the z-th rightmost boundary point and Qi be the position of the z-th rightmost 
generator (z = 1, . . . , 5). It is easy to verify that the only admissible configuration for the boundary points in order 



March 30, 2009 



DRAFT 



ffiEE TRANSACTIONS ON AUTOMATIC CONTROL, SUBMITTED FOR PUBLICATION AS A FULL PAPER. 



17 




Fig. 4. Example of non-existence of an equitable Voronoi diagram on a line. The above tessellation is an equitable partition, but not a Voronoi 
diagram. 

to obtain an equitable Voronoi diagram is the one depicted in Fig. [4] Now, by the Perpendicular Bisector Property, 
we require: 

{g3-b2 = b2- 92 

Therefore, we would require 94—92 = 2(63 — 62) = 1-2; this is impossible, since 92 G [0.1, 0.2] and 94 G [0.8, 0.9]. 

VI. Distributed Algorithms for Equitable Partitions with Special Properties 

The gradient descent law ([17]), although effective in providing convex equitable partitions, can yield long and 
"skinny" subregions. In this section, we provide spatially-distributed algorithms to obtain convex equitable partitions 
with special properties. The key idea is that, to obtain an equitable power diagram, changing the weights, while 
maintaining the generators fixed, is sufficient. Thus, we can use the degrees of freedom given by the positions of 
the generators to optimize secondary cost functionals. Specifically, we now assume that both generators' weights 
and their positions obey a first order dynamical behavior 

Wi = uf ^ 
Define the set 

S= e {A xR)"^ \9i 9j for Silli^ j.SindXv, >0 Vi G /^n}. (18) 

The primary objective is to achieve a convex equitable partition and is captured, similarly as before, by the cost 
function Hy : S \-^Ryo 

m 

Hv{Gw) = Xl^v.W)- 

We have the following 
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Theorem 6.1: Given a measure A, the function Hy is continuously differentiable on 5, where for each i G 
{l,...,m} 

dH 1 1 1 ^^^^ 

Furthermore, the critical points of Hy are generators' positions and weights with the property that all power cells 
have measure equal to A^/m. 

Proof: The proof of this theorem is very similar to the proof of Theorem |5.6[ we omit it in the interest of 
brevity (the derivation of the partial derivative can be found in ifTTl ). ■ 



Notice that the computation of the gradient in Theorem 6.1 is spatially distributed over the power-Delaunay graph. 
For short, we define the vectors v±qh. = i^^- 

Three possible secondary objectives are discussed in the remainder of this section. 

A. Obtaining Power Diagrams Similar to Centroidal Power Diagrams 
Define the mass and the centroid of the power cell F^, z G /rn. as 

My. — I \(x) dx, Cy = — — / xX(x) dx. 
' Jv, ' ' My, Jy^ ^ 

In this section we are interested in the situation where gi = Cy., for all i e Im- We call such a power diagram a 
centroidal power diagram. The main motivation to study centroidal power diagrams is that, as it will be discussed 



in Section VI-C the corresponding cells, under certain conditions, are similar in shape to regular polygons. 

A natural candidate control law to try to obtain a centroidal and equitable power diagram (or at least a good 
approximation of it) is to let the positions of the generators move toward the centroids of the corresponding regions 
of dominance, when this motion does not increase the disagreement between the measures of the cells (i.e., it does 
not make the time derivative of Hy positive). 

First we introduce a saturation function as follows: 

' for X < , 

exp(^-p^^ forx>0, G M>o. 
Moreover, we define the vector vc^g^ = Cy. — gi. Then, we set up the following control law defined over the set 
S, where we assume that the partition V{Gw) = {^i? • • • , Kn} is continuously updated, 

dHy 



e{x) = 



OWi 

2 

gi = — arctan 

TT 



a 



(20) 



In other words, gi moves toward the centroid of its cell if and only if this motion is compatible with the 
minimization of Hy. In particular, the term arctan ^ || is needed to make the right hand side of ( |2Q| ) 
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C^, while the term &{vc,gi '^-oh ) needed to make the right hand side of ( [2Q| compatible with the minimization 
of Hv- To prove that the vector field is it is simply sufficient to observe that it is the composition and product 
of functions. Furthermore, the compatibility condition of the flow ( [2Q| ) with the minimization of Hy stems from 
the fact that = as long as vc,gi ' "^-qh — ^' presence of 6(-)- Notice that the computation of the 

right hand side of ( [20| ) is spatially distributed over the power-Delaunay graph. 

As in many algorithms that involve the update of generators of Voronoi diagrams, it is possible that under control 
law ( [2Q| there exists a time t* and j G Im such that gi{t*) = gj{t*). In such a case, either the power diagram is 
not defined (when Wi{t*) = Wj{t^)), or there is an empty cell (Wi{t*) / Wj{t^)), and there is no obvious way to 
specify the behavior of control law ( [2Q| for these singularity points. Then, to make the set S positively invariant, 
we have to slightly modify the update equation for the positions of the generators. The idea is to stop the positions 
of two generators when they are close and on a collision course. 

Define, for A e M>o, the set 

M,(G, A) = {g, eG\ \\g, - g^ \\ < A, g, + g,}. 

In other words. Mi is the set of generators' positions within an (Euclidean) distance A from gi. For b G M>o, 
^ < A, define the gain function i\){p, i3) : [0, A] x [0, 27r] ^ R>o (see Fig. [S]): 

^ if (5 < p < A and < I? < TT, 

^(1 + sini9) - sin^ if S < p < A and tt < ^ < 27r, 
if p<S and < 1? < TT, 

-fsin^^ if p<S and7r<^<27r, 

It is easy to see that ?/^(-) is a continuous function on [0, A] x [0, 27r] and it is globally Lipschitz there. Function 
?/;(•) has the following motivation. Let p be equal to — and let Vx be a vector such that the tern {vx,{gj — 
gi)^'^xX {gj —gi)} an orthogonal basis of R^, co-orientied with the standard basis. In the Figure [sj Vx corresponds 
to the X axis and gj — gi corresponds to the y axis. Then the angle is the angle between Vx and vcg^, where 
< i^ij < 27r. If p < ^ and < i^ij < tt, then gi is close to gj and it is on a collision course, thus we set the 
gain to zero. Similar considerations hold for the other three cases; for example, if p < ^ and tt < i^ij < 27r, the 
generators are close, but they are not on a collision course, thus the gain is positive. 
Thus, we modify control law ( [20| ) as follows: 

- — ,/Cent,w 

OWi 

(22) 

If Mi(G, A) is the empty set, then we have an empty product, whose numerical value is 1. Notice that the right 
hand side of ([22]) is Lipschitz continuous, since it is a product of functions and Lipschitz continuous functions, 
and it can be still computed in a spatially-distributed way (in fact, it only depends on generators that are neighbors 
in the power diagram, and whose positions are within a distance A). One can prove the following result. 



2 

— arctan 

TT 



\\v. 



a 
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Zero gain 



Fig. 5. Gain function used to avoid that the positions of two power generators can coincide. 



Theorem 6.2: Consider the vector field on S defined by equation ( [22] ). Then generators' positions and weights 
starting at t = at Gw{^) ^ and evolving under ([22]) remain in S and converge asymptotically to the set of 
critical points of the primary objective function Hy (i.e., to the set of vectors of generators' positions and weights 
that yield an equitable power diagram). 



Proof: The proof is virtually identical to the one of Theorem |5.9[ and we omit it in the interest of brevity. We 
only notice that Hy is non-increasing along system trajectories 



^dHy . 



dHy 

dwi 




fdHv 

\ dwj 



< 0. 



Moreover, the components of the vector field ([22]) for the position of each generator are either zero or point toward 
A (since the centroid of a cell must be within A); therefore, each generator will remain within the compact set A. 



B. Obtaining Power Diagrams ''Close'' to Voronoi Diagrams 

In some applications it could be preferable to have power diagrams as close as possible to Voronoi diagrams. 
This issue is of particular interest for the setting with non-uniform density, when an equitable Voronoi diagram 
could fail to exist. The objective of obtaining a power diagram close to a Voronoi diagram can be translated in the 
minimization of the function K : ^ R>o: 



m 



when Wi = for all i e Im, have K = and the corresponding power diagram coincides with a Voronoi 
diagram. To include the minimization of the secondary objective K, it is natural to consider, instead of ( [TT] ), the 
following update law for the weights: 

dHy dK dHy 



Wi 



dwi 



dwi 



dwi 



Wi. 



(23) 
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However, Hy is no longer a valid Lyapunov function for system ( [23] ). The idea, then, is to let the positions of the 
generators move so that • gi — = 0. In other words, the dynamics of generators' positions is used to 

compensate the effect of the term —Wi (present in the weights' dynamics) on the time derivative of Hy- 
Thus, we set up the following control law, with ^i, £2 and £3 positive small constants, £2 > £1, 

dHv 



dwi 



iS^ts,,s2 [W'^dHi 11) sato,£3 (dist(^i, dVi) 



dH V ~ ^^^^ 

The gain sat£^^£2 (iba^Jl) needed to make the right hand side of p4| ) Lipschitz continuous, while the gain 
sato,£3 (^dist{g i^dVi)^ avoids that generators' positions can leave the workspace. Notice that the computation of the 
right hand side of p4| ) is spatially distributed over the power-Delaunay graph. 

As before, it is possible that under control law p4| ) there exists a time t* and j G Im such that gi{t*) = gj(t*). 
Thus, similarly as before, we modify the update equations ([24]) as follows 

= - ^isat^,,£2 (iba^Jl) sato,£3 (dist(^i,51/i)) • W ^lj{\\gj - gi\\, ^ij^ = ^if''"^, 

OH ^ ^ 



where is defined as in Section 



VI-A 



with Wi^^VQfj^ replacing vc,g,. 
One can prove the following result. 

Theorem 6.3: Consider the vector field on S defined by equation ( [25] ). Then generators' positions and weights 
starting at t = at Gw{^) ^ and evolving under ( [25] ) remain in S and converge asymptotically to the set of 
critical points of the primary objective function Hy (i.e., to the set of vectors of generators' positions and weights 
that yield an equitable power diagram). 

Proof: Consider Hy as a Lyapunov function candidate. First, we prove that set S is positively invariant with 
respect to ( [25] ). Indeed, by definition of ( [25] ), we have g^ 7^ gj for i / j for alH > (therefore, the power diagram 
is always well defined). Moreover, it is straightforward to see that Hy < 0. Therefore, it holds 

K^G^it)) < Hy{Gw{t)) < Hy{Gwm. i ^ Im. t > 0. 

Since the measures of the power cells depend continuously on generators' positions and weights, we conclude that 
the measures of all power cells will be bounded away from zero. Finally, since = on the boundary of A for all 
i e Im, each generator will remain within the compact set A. Thus, generators' positions and weights will belong 
to S for all t > 0, that is, Gw{t) G 5' Vt > 0. 

Second, as shown before, Hy : S ]R>o is non-increasing along system trajectories, i.e., IIy{Gw) < in 
Third, all trajectories with initial conditions in S are bounded. Indeed, we have already shown that each generator 
remains within the compact set A under control law ([25]). As far as the weights are concerned, we start by noticing 
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that the time derivative of the sum of the weights is 

i=l gjeMi{G,A) 



since, similarly as in the proof of Theorem 5.9 XI = 0. Moreover, the magnitude of the difference between 



any two weights is bounded by a constant, that is, 

I'^^i — Wj\ ^ <^ for all j e Im'', (26) 

if, by contradiction, the magnitude of the difference between any two weights could become arbitrarily large, the 
measure of at least one power cell would vanish, since the positions of the generators are confined within A. 
Assume, by the sake of contradiction, that weights' trajectories are unbounded. This means that 

Vi? > 3t > and 3 j G Im such that \wj{t)\ > R. 

For simplicity, assume that Wi{0) = for all i e Im (the extension to arbitrary initial conditions in S is 
straightforward). Choose R = 2ma and let t2 be the time instant such that \wj{t2) \ = R, for some j ^ Im- Without 
loss of generality, assume that Wj{t2) > 0. Because of constraint ([26]), we have Yl^i ^1(^2) > f m(3m + 1). Let 
ti be the last time before ^2 such that Wj{t) = ma; because of continuity of trajectories, ti is well-defined. Then, 
because of constraint ([26]), we have (i) X^I^i '^i(^i) ^ f m(3m — 1) < X^I^i '^^(^2), <^nd (ii) ^^'^^ ^''^^^ ^ for 
t e [ti, ^2] (since Wj{t) > ma for all t e [ti, ^2] and Eq. ( [26] ) imply min^^/^ Wi{t) > for all t e [ti, ^2]); thus, 
we get a contradiction. 



Finally, by Theorem 6.1 Hy is continuously differentiable in S. Hence, by invoking the LaSalle invariance 
principle (see, for instance, |6|), under the descent flow ( [25] ), generators' positions and weights will converge 
asymptotically to the set of critical points of Hy, that is not empty by Theorem 



5.1 



C. Obtaining Cells Similar to Regular Polygons 

In many applications, it is preferable to avoid long and thin subregions. For example, in applications where a 
mobile agent has to service demands distributed in its own subregion, the maximum travel distance is minimized 
when the subregion is a circle. Thus, it is of interest to have subregions whose shapes show circular symmetry, i.e., 
that are similar to regular polygons. 

Define, now, the distortion function Ly ' {A x R)^ \ Fcoinc ^ ^>o' Yl^i Jy Ik ~ gi\\'^X{x)dx (where Vi is 
the i-th cell in the corresponding power diagram). In [|18i it is shown that, when m is large, for the centroidal 
Voronoi diagram (i.e., centroidal power diagram with equal weights) that minimizes Ly, all cells are approximately 



congruent to a regular hexagon, i.e., to a polygon with considerable circular symmetry (see Section VII for a more 
in-depth discussion about circular symmetry). 

Indeed, it is possible to obtain a power diagram that is close to a centroidal Voronoi diagram by combining 
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control laws ([22]) and ( [25] ). In particular, we set up the following (spatially-distributed) control law: 

or,w 

(27) 



cent,g 



Combining the results of Theorem 6.2 and Theorem 6.3 we argue that with control law ( [27] ) it is possible to obtain 



equitable power diagrams with cells similar to regular polygons, i.e. that show circular symmetry. 

VII. Simulations and Discussion 

In this section we verify through simulation the effectiveness of the optimization for the secondary objectives. Due 
to space constraints, we discuss only control law ( [27] ). We introduce two criteria to judge, respectively, closeness 
to a Voronoi diagram, and circular symmetry of a partition. 

A. Closeness to Voronoi Diagrams 

In a Voronoi diagram, the intersection between the bisector of two neighboring generators gi and gj, and the 
line segment joining gi and gj is the midpoint ^JJ^ = {gi -\- gj)/2. Then, if we define g^j"^ as the intersection, in 
a power diagram, between the bisector of two neighboring generators {gi^Wi) and {gj^Wj) and the line segment 
joining their positions gi and gj, a possible way to measure the distance of a power diagram from a Voronoi 
diagram is the following: 

"^-^1^1^ 0.57., ' ^ ^ 

where N is the number of neighboring relationships and, as before, jij = \\gj — gi\\. Clearly, if the power diagram 
is also a Voronoi diagram (i.e., if all weights are equal), we have r] = 0. We will also refer to r] as the Voronoi 
defect of a power diagram. 

B. Circular Symmetry of a Partition 

A quantitative manifestation of circular symmetry is the well-known isoperimetric inequality which states that 
among all planar objects of a given perimeter, the circle encloses the largest area. More precisely, given a planar 
region V with perimeter py and area |V|, then — 47r|V| > 0, and equality holds if and only if F is a circle. 
Then, we can define the isoperimetric ratio as follows: Qv = ^^^^^ ; by the isoperimetric inequality, Qy < 1, 

Pv 

with equality only for circles. Interestingly, for a regular n-gon the isoperimetric ratio Qn is Qn = ntln^ ' which 
converges to 1 for n ^ oo. Accordingly, given a partition A = {A^}^^, we define, as a measure for the circular 
symmetry of the partition, the isoperimetric ratio Qj[ = QAi • 

C. Simulation Results 

We simulate ten agents providing service in the unit square A. Agents' initial positions are independently and 
uniformly distributed over A, and all weights are initialized to zero. Time is discretized with a step dt = 0.01, 
and each simulation run consists of 800 iterations (thus, the final time is T = 8). Define the area error e as the 
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TABLE I 

Performance of control law pT) . 



A 


E [e] max e 


lEW 


max?7 




minQv 


unif 


3.810-4 0.016 


0.01 


0.03 


0.73 


0.66 


gauss 


310-3 5.310-3 


0.02 


0.04 


0.75 


0.69 



difference, at T = 8, between the measure of the region of dominance with maximum measure and the measure of 
the region of dominance with minimum measure. 

First, we consider a measure A uniform over A, in particular A = 1. Therefore, we have A^ = 1 and agents 
should reach a partition in which each region of dominance has measure equal to 0.1. For this case, we run 50 
simulations. 

Then, we consider a measure A that follows a gaussian distribution, namely \{x^y) = e-^((^-^-^)^+(^-^-^)^), 
(x, y) G A, whose peak is at the north-east corner of the unit square. Therefore, we have A^ ^ 0.336, and vehicles 
should reach a partition in which each region of dominance has measure equal to 0.0336. For this case, we run 20 
simulations. 

Table |l| summarizes simulation results for the uniform A (A=unif) and the gaussian A (A=gauss). Expectation and 
worst case values of area error e, Voronoi defect r] and isoperimetric ratio Qy are with respect to 50 runs for uniform 
A, and 20 runs for gaussian A. Notice that for both measures, after 800 iterations, (i) the worst case area error is 
within 16% from the desired measure of dominance regions, (ii) the worst case r] is very close to 0, and, finally, 
(iii) cells have, approximately, the circular symmetry of squares (since Q4 ^ 0.78). Therefore, convergence to a 
convex equitable partition with the desired properties (i.e., closeness to Voronoi diagrams and circular symmetry) 
seems to be robust. Figure |6] shows the typical equitable partitions that are achieved with control law ( [27] ) with 10 
agents. 

VIII. Application and Conclusion 
In this last section, we present an application of our algorithms and we draw our conclusions. 

A. Application 

A possible application of our algorithms is in the Dynamic Traveling Repairman Problem (DTRP). In the DTRP, 
m agents operating in a workspace A must service demands whose time of arrival, location and on-site service 
are stochastic; the objective is to find a policy to service demands over an infinite horizon that minimizes the 
expected system time (wait plus service) of the demands. There are many practical settings in which such problem 
arises. Any distribution system which receives orders in real time and makes deliveries based on these orders (e.g., 
courier services) is a clear candidate. Equitable partitioning policies (with respect to a suitable measure A related 
to the probability distribution of demand locations) are, indeed, optimal for the DTRP when the arrival rate for 
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(a) Typical equitable partition of A for A(a::, y) = 1. (b) Typical equitable partition of A for X(x,y) = 

g-5((cc-0.8)2 + (y-0.8)2)^ 

Fig. 6. Typical equitable partitions achieved by using control law pT) . The yellow squares represent the position of the generators, while the 
blue circles represent the centroids. Notice how each bisector intersects the line segment joining the two corresponding power neighbors almost 
at the midpoint; hence both partitions are very close to Voronoi partitions. Compare with Fig. [2] 



service demands is large enough (see Q, (191, COl)- Therefore, it is of interest to combine the optimal equitable 
partitioning policies in |19 | with the spatially-distributed algorithms presented in this paper. 

The first step is to associate to each agent i a virtual power generator (virtual generator for short) (gi^Wi). We 
define the region of dominance for agent i as the power cell Vi = Vi{Gw)^ where Gw = (^(^i, ^i), • • • , (^m, ^m)^ 
(see Fig. |7(a) ). We refer to the partition into regions of dominance induced by the set of virtual generators Gw as 
V{Gw)' A virtual generator (gi^Wi) is simply an artificial variable locally controlled by the i-th agent; in particular, 
gi is a virtual point and Wi is its weight. Virtual generators allow us to decouple the problem of achieving an equitable 
partition into regions of dominance from that of positioning an agent inside its own region of dominance. 

Then, each agent applies to its virtual generator one of the previous algorithms, while it performs inside its region 



of dominance the optimal single-agent policy described in |1| (see Fig. |7(b)| ). 

Notice that, since each agent is required to travel inside its own region of dominance, this scheme is inherently 
safe against collisions. 



B. Conclusion 

We have presented provably correct, spatially-distributed control policies that allow a team of agents to achieve 
a convex and equitable partition of a convex workspace. We also considered the issue of achieving convex and 
equitable partitions with special properties (e.g., with hexagon-like cells). Our algorithms could find applications in 
many problems, including dynamic vehicle routing, and wireless networks. This paper leaves numerous important 
extensions open for further research. First, all the algorithms that we proposed are synchronous: we plan to devise 
algorithms that are amenable to asynchronous implementation. Second, we envision considering the setting of 
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(a) Agents, virtual generators and regions of dominance. (b) Each agent services outstanding demands inside its 

own region of dominance. 

Fig. 7. Spatially-distributed algorithms for the DTRP. 



Structured environments (ranging from simple nonconvex polygons to more realistic ground environments). Finally, 
to assess the closed-loop robustness and the feasibility of our algorithms, we plan to implement them on a network 
of unmanned aerial vehicles. 
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Appendix 

Proof of Theorem 5.16- The proof mainly relies on II2TII . Let v be the unit vector considered in the definition 



of the Unimodal Property. Then, there exist unique values sq < si < ■ ■ ■ < Sm such that sq = mf{s;A^ 7^ 0}, 
Sm = sup{s; 0}, and 

k 

>^{xeA;vx<Sk} = —^A, /c = 1, . . . ,m - 1. (29) 

Consider the intervals li = [si-i,Si], i e Im- We claim that one can choose points gi = Uv G R^, i ^ Im such 
that ti e li and the corresponding Voronoi diagram is 

Ai = {x e A; \\x - gi\\ min \\x - gk\\} 

^ (30) 
= {x e A; V -x e [si-i,Si]}. 

Together, Eq. ([29]) and Eq. ( [3Q| ) yield the desired result. 

Since, by assumption, A enjoys the Unimodal Property, there exists an index G {1, . . . , m} such that the length 

of the intervals 1^ = decreases as i ranges from 1 to hz, then increases as i ranges from tz to m. Let 

Ik = [sk-i^s^] be the smallest of these intervals, and define 

By induction, for i increasing from to m — 1, define t^+i be the symmetric to ti with respect to 5^, so that 

tij^i = 2si — ti i = K, -\- 1, . . . ,m — 1. 
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Since the length of is larger than the length of we have 

U e li ^U^i e /i+i. (31) 

Similarly, for i decreasing from to 1, we define 

ti—i 25^—1 ^i, i = hv^ tv 1, . . . , 2. 
Since the interval li-i is now larger than the interval we have 

U e li ^ U-i e li-i. (32) 

By Eqs. ([3T])-([32|) imply ti e li for alH = 1, . . . , m. Hence the second equality in Eq. ( [3Q| ) holds. ■ 
We now specialize the theorem to the case when A is convex. 

Corollary 8.1: Let A C be a compact, convex set, and A be constant on A. Then for every m > 1 there exist 
points gir ■ ■ ^Qm all in the interior of A, such that the corresponding Voronoi diagram is equitable. 

Proof: Notice that every compact convex set enjoys the Unimodal Property, with an arbitrary choice of the unit 
vector V. By compactness, there exist points b e A such that ||^ — a|| = max^^^'GA Ik ~ ^11- a translation of 
coordinates, we can assume a = 0. Choose v = Then the previous construction yields an equitable Voronoi 

diagram generated by m points gi = tiV all in the interior of A. ■ 



Proof of Lemma 5.14' By Theorem 5.9 and by its very definition W*{G) is the zero of the vector field 

V 

dwi 



^^^-{W{t)). Now let us denote with 



K{W,G)=-^{W), 

the corresponding continuous function, viewed as a function of two independent set of variables, namely the 
weights {wi^ . . . ^Wn) = W and the non-degenerate vector of generators' locations G. In order to prove that 



the assignment G W*{G) is continuous, notice that by Theorem 5.6 the function K{W^G) is identically 
zero when restricted to the graph of W*, namely K{W*{G)^G) = 0. The function W* is continuous iff it is 
continuous in each of its argument. Fix, first, a generator ^ dT and consider for any v G M^, the variation 
{gii " ' igi-i^gi + hv, gi-^i^ . . . , gm)- Since g^ ^ dV, there always exists an e > 0, depending on gi and v, 
such that for any h with < h < e, {gi, . . ..gi^i.gi + hv.gi^i, . . .,gm) ^ F. Now K{W*{gi, . . .,gi-i,gi + 
hv, ^i+i, . . . , gm), (^1, . . . , gi-i,gi + hv, gi^i, . . . , gm)) = for any < /i< e by definition. Therefore, taking the 
limit for ^ 0+, we still get zero. On the other hand, since K is continuous, we can take the limit inside K and 
we get 

K{ lim W*(^i,...,^i_i,^i + /iv,^i+i,...),(^i,...,^i_i,^i,^i+i,...)) = 0. 

Therefore, we have that lim/,^o+ • • • , gi-i,gi^hv, ^^+1, . . . , is equal to W*(^i, . . . , ^i+i, . . . , 

by the uniqueness in ft of the value of for which, given G, the function K vanishes. ■ 
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